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Abstract 

A self-consistent equation to derive a discreteness-induced stochastic steady state is presented for reaction- 
diffusion systems. For this formalism, we use the so-called Kuramoto length, a typical distance over which 
a molecule diffuses in its lifetime, as was originally introduced to determine if local fluctuations influence 
globally the whole system. We show that this Kuramoto length is also relevant to determine whether the 
discreteness of molecules is significant or not. If the number of molecules of a certain species within the 
Kuramoto length is small and discrete, localization of some other chemicals is brought about, which can 
accelerate certain reactions. When this acceleration influences the concentration of the original molecule 
species, it is shown that a novel, stochastic steady state is induced that does not appear in the continuum 
limit. A theory to obtain and characterize this state is introduced, based on the self-consistent equation for 
chemical concentrations. This stochastic steady state is confirmed by numerical simulations on a certain 
reaction model, which agrees well with the theoretical estimation. Formation and coexistence of domains 
with different stochastic states are also reported, which is maintained by the discreteness. Relevance of our 
result to intracellular reactions is briefly discussed. 

1 Introduction 

Chemical reaction dynamics are often studied with the use of rate equations for chemical concentrations. For 
this approach, the number of molecules is assumed to be large, which validates the continuum description. 
However, in a biological system such as a cell, the number of molecules within is sometimes rather small. Then 
the validity of continuum description by the rate equations is not evident. This problem of smallness in molecule 
number is not restricted in biology. Following recent advances in nanotechnology, reactions in a micro-reactor 
are studied experimentally, where the number of molecules in concern is quite small. This is also true in some 
surface reaction of absorbed chemicals. 

Here we are interested in the effect of such smallness in molecule number. Of course, one straightforward 
consequence of the smallness in the number is the large fluctuations in the concentration. Indeed, the fluctuations 
around the continuous rate equation can be discussed by stochastic differential equation ^ |2] . State change by 
noise has been studied as noise-induced transitions 0, noise- induced order py, stochastic resonance [2], and so 
forth. The use of stochastic differential equation, as well as its consequence, has been investigated thoroughly. 
If the number of molecules is much smaller and can reach 0, however, another effect of "smallness" is expected, 
that is the discreteness in the number. Our concern in the present paper is a drastic effect induced by such 
discreteness in the molecule number. 

Previously we have discovered a transition of a chemical state, induced by the discreteness in the molecule 
number, i.e., the effect of the number of molecules 0, 1, • • • joj. The transition is termed as discreteness- induced 
transition (DIT). This transition is not explained by the stochastic differential equation approach. Rather, 
discreteness, in conjunction with the stochastic effect, is essential. In particular, we studied a system with 
autocatalytic reaction in a well stirred container. As the volume of the container decreases and the molecule 
number decreases, a transition to a novel state with symmetry breaking occurs, that does not appear either 
in the continuous rate equation or in its Langevin version. Here the transition occurs, when the number of 
molecule flow from environment to the reactor is discrete, in the sense that it is less than one on the average, 
within the average reaction time. Indeed, to the discreteness-induced transition, relevant is not the molecule 
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number itself but the discreteness in the number of some molecular process (e.g., flow of molecule into the 
system) within the average time scale of some other reaction process. 

On the other hand, in a spatially extended system with reaction and diffusion, the total number of molecules 
(and molecular events) increases with the system size, and is not small. Instead, the number of molecules (or 
events), not in the total system but within the size of an "effective length", is relevant to determine the 
discreteness effect. Then, we need to answer what this effective length is. In 0, we have proposed that the 
so-called Kuramoto length gives an answer to it. 

Kuramoto length £k is defined as the average length that a molecule diffuses within its lifetime, i.e., before 
it makes reaction with other molecules J] [5] [5| . In the seminal papers [5] |5] , Kuramoto has shown that whether 
the total system size is larger than this length or not provides a condition to guarantee the use of the reaction- 
diffusion equation. When the system size (length) is smaller than local fluctuations rapidly spread over the 
system. Contrastingly, if the system size is much larger than £k , distant regions fluctuate independently, and 
the system is described by local reaction process and diffusion, validating the use of reaction-diffusion equation. 

For example, consider the reaction 

A X, 2X B. 

If the concentration of chemical A is set to be constant, the chemical X is produced at the constant rate k, 
while it decays with the reaction 2X — > B at the rate k' . The average concentration of X at the steady state is 
(X) = y/kA/2k', where A is the concentration of the chemical A. Thus the average lifetime of X at the steady 
state is estimated to be r = l/(2k'(X)) = l/\/2kk'A. If X molecules diffuse at the diffusion constant D in 
one-dimensional space, the typical length over which an X molecule diffuses in its lifetime is estimated to be 

£ K = V2#7, (1) 

which gives the Kuramoto length. 

In these works, it is assumed that the average distance between molecules is much smaller than £k, and 
there is a large number of molecules within the region of the length Ik- Thus the concentration of the chemical 
X can be regarded as a continuous variable. Hence the continuum description is valid. However, if the average 
distance between molecules is comparable to or larger than £k, local discreteness of molecules may not be 
negligible. 

For example, consider a chemical species Xj, whose Kuramoto length is given by £j. Then we consider 
discreteness of molecule species Xj that produces this chemical Xj, i.e., the case that average number of Xj 
is less than 1 within the area of the Kuramoto length £j. With this setting, molecules Xj, produced by Xj 
molecules, will be localized around them, as the average distance between Xj molecules is larger than the 
Kuramoto length of Xj. Then, this localization of the chemical may drastically alter the total rate of the 
reactions, if reactions with 2nd or higher order of Xj are involved, as will be shown later. In the present paper, 
following we pursue the possibility that discreteness of some molecules within Kuramoto length of some 
other molecules may drastically change the steady state of the system, as in DIT previously studied. 

In section 2, wc discuss a general condition for the amplification of some reaction by such discreteness. Then 
by introducing a self-consistent equation for the rate of this amplification, we demonstrate the existence of stable 
stochastic steady state (SSS), that never appears in the continuum description. In section 3, we numerically 
study a specific chemical reaction model with three components, to show the validity of this self-consistent 
theory for SSS. In section 4, domain formation with this SSS is presented, as a novel possibility for pattern 
formation in reaction-diffusion system. Discussion is given in section 5, with possible applications to biological 
problems. 

2 Steady state induced by discreteness of molecule, with amplifica- 
tion of some reaction: self-consistent analysis 

Consider a reaction system consisting of several molecule species X m (to = with chemical reaction 

and diffusion. The system can involve catalytic reactions of higher-order catalysis or autocatalysis. Some other 
molecules (e.g., resource chemicals) are supplied externally, involved in the reaction among X m , so that the 
noncquilibrium condition is sustained. So far the system in concern is rather general chemical reaction system 
with diffusion. 
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Now, we take a pair of molecule species, Xi and Xj, where Xj is produced by Xi, and study how the 
discreteness of the molecule Xi can alter the steady state from the continuum limit case. To discuss the 
discreteness effect, wc consider the case that the molecule Xj is localized around the molecule Xi. (Recall the 
molecule Xj is produced by Xi.) In order for this localization to work, the average length that the molecule 
Xj travels within its lifetime should be smaller than the average distance of the molecules X$. In other words, 
the average number of the molecules Xi within the domain of the Kuramoto length £j is less than 1. (see Fig. 
1). Here, the lifetime of the molecule Xj is determined by the collision with some other molecule species whose 
density is not low. Hence, the Kuramoto length of the Xj molecule, determined as in §1, is given independently 
of the concentrations of the molecules Xi and Xj. 

Now, to alter drastically the steady state by discreteness, the localization of Xj molecule has to change 
concentrations of some other molecules, as compared with the case of homogeneous distribution of Xj. This is 
possible if there is a higher-order reaction such as mXj + X q — > X p , because the probability of such reaction is 
amplified by localization of Xj molecules in space. To compute this acceleration, we calculate the average of cj 1 , 
where Cj is the concentration of Xj , and compute the degree of amplification a from that for the homogeneously 
distributed case. In the calculation, we assume that Xj is localized around the Xi molecules, with a width of 
lj 1 the Kuramoto Length of Xj (which is shorter than the average distance between X{ molecules). 

Assuming that the distribution of Xj is represented by the continuous concentration Cj (x), a can be expressed 

as 

where V is the size (volume) of the system. 

For simplicity, we assume that Xi is randomly distributed over d-dimensional space, and the distribution of 
Xj is given by a d-dimensional Gaussian distribution with a standard deviation £j around each Xi molecule, 
such as 

where x~k is the position of each Xi molecule. Now Cj(x) is the sum of Pk(x); thus, (cj) — Ci since J pkdx = 1. 
For the case with m = 2 and sufficiently large V, 



(s 2 ) = {(Ep k ) 2 ) 



since Xi is randomly distributed. With eq. J5J, 

a = ({Cj) 2 + (2^ij)- d {cj))/{cj) 2 

= 1 + {2^ij)- d { Cj )- x = 1 + (20F)- W 



Thus, we obtain the acceleration factor 



a = 1 + ~, — 7=7-, 1- (3) 

{2^y art v ; 



In the same manner, for m = 3, 



and generally, for Cjff <C 1, 



a = 1 + =T3 7 H 7= ; — . (4) 

{2^Y (2^y (^) 2 



a^m-^^f^cdf) m ■ (5) 

As shown, the reaction can be drastically amplified as the number of Xi molecules within a volume of the 
Kuramoto length (cilj) is much smaller than 1. 

So far we have shown that for reaction system involving the process from Xi to Xj, the discreteness can 
alter the concentration of some chemicals drastically, if (1) the density of Xi molecule is so low that the number 
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is discrete within the size of the Kuramoto length of Xj molecule and (2) there is a high order (higher than 
linear) reaction with regards to Xj. 

Next, to confirm that this acceleration of reaction alters the steady state from the continuum case, we need 
to check if the condition for the discreteness is sustained under the above amplification of concentration of some 
chemicals, as a steady state solution. Hence we study some feedback from the concentration of X p to Xi that 
is generated by some reaction path(s). If Xi is produced or catalyzed by X p , the concentration of Xi depends 
on that of X p , c p . With such feedback, the change of concentration Cj is given by some function F(cj, c p ), while 
the change of the concentration Cj depends on Cj, and is given by some function as G(ci). Since c p is a function 
of a(ci), F(cj, c p ) is rewritten as F(cj, a(cj)). Hence the concentrations of Xi and Xj molecules must satisfy 

dci/dt — F(cj,a(ci)); dcj/dt = G(ci) (6) 

(F and G may have dependence on other concentrations or reaction rates, e.g., G{ci) can also depend on Cj or 
a). The steady state solution is obtained by setting the right hand of these equations as 0. 

Note that the solution with a = 1 corresponds to the continuum case, given by the standard rate equation. 
For some case, this is the only solution for the concentrations. For some other cases, however, there is some 
other solution(s) with a > 1. This is a solution with the amplification by localization of molecules due to 
the discreteness of the Xi molecule. If the concentration of Cj obtained from this solution satisfies Citj < 1, 
this discreteness-induced solution is self-consistent. Furthermore, the stability of this solution is computed by 
linearizing the solution around this fixed point. If this solution is linearly stable, stability of this novel steady 
solution is assured, which does not exist in the continuum description (or in its Langevin equation version). 
We call the state represented by this solution as stochastic steady state (SSS), as it is sustained stochastically 
through discreteness in molecule numbers. We will show an explicit example of this SSS in the next section. 

Self-consistent solution involving the change of Kuramoto-length 

So far we have assumed that the Kuramoto-length tj of the Xj molecule is constant. This is true as long 
as the concentration of the molecule relevant to the decomposition or transformation of the Xj molecule is 
constant. However, if the concentration of the chemical that is relevant to the determination of tj depends on 
the concentration of either Xi, Xj, or X p , the Kuramoto length, as well as a, depends on it. Accordingly, in 
eq. ©, we need to regard tj in a as a variable that depends on cither a, Cj or a. With the inclusion of the 
dependence, we again obtain a self-consistent solution, to get the concentrations Cj and Cj (and accordingly a 
and tj). If there is a stable solution with a > 1 and Citj < 1, then we get a SSS as a self-consistent solution 
both on a and tj. We will discuss a related example in §4, where two solutions with a = 1 and a > 1 coexist 
in space, and form a domain structure. 

Combination of several processes 

So far we have discussed a simple case of discreteness-induced state. The discussion with the use of amplification 
factor a, however, is generalized to include temporally or spatially dependent solutions of c$ and Cj, with 
temporal (or spatial) dependence of a. This solution represents an average behavior longer than the time scale 
for stochastic collisions or longer scale than tj. With this extension, we can discuss discreteness- induced rhythm 
or pattern, that is stochastically sustained. 

Such spatiotemporal dynamics can often appear in a reaction network of several molecules, with two or 
more pairs of discreteness in number. For example, consider reactions Xn — > {produces) — > Xji, Xji — ► 
(produces with high order reaction) — > X p \, X p \ — > (produces) — > X^; and Xi2 — » (produces) — > Xji, 
Xji — > (produces with high order reaction) —> X p i, X p i — > (produces) — » Xn (see Fig. EJ, where we assume 
that the density en, en, of Xn and Xn molecules are so low that Cnlj X < 1 and cntj 2 < 1, respectively for the 
Kuramoto lengths of Xji and Xji. Then, following the scheme we discussed, we get a coupled equation for the 
concentrations of en, Cj\, en, Cji with two amplification factors a,\ and cti. In general, there may be a time- 
or space- dependent solution (by including diffusion term with much longer spatial scale), that leads to a novel 
stochastic pattern or rhythm. Explicit examples for such case will be discussed in future. 
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(5) 



Figure 1: Schematic representation for the mechanism to produce SSS. Assume that (1) the chemical Xj is 
produced by Xi, and (2) lj, the Kuramoto length of Xj, is shorter than Aj(= c^ 1 ), the average distance of Xi 
molecules (i.e., Ci£j < 1). Then, Xj is localized around Xi molecules. (3) If there is m-th order reaction of Xj 
(to > 1), the rate of the reaction is proportional to c™ 1 ; hence, the reaction is accelerated. (4) Additionally, if 
the reaction promotes (directly or indirectly) the production of Xi, the acceleration of the reaction may cause 
increase of Cj. (5) On the other hand, if Cj is high (c^ « 1 or > 1), Xj is almost uniformly distributed; thus, 
the acceleration is weak, and the production of X, is degraded. Hence, there is a stable steady state of a at an 
intermediate value. 
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Figure 2: Schematic diagram of an example of reaction cascade. 
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3 Specific example of Stochastic Steady State 



To confirm our theoretical estimation for SSS, we have adopted a simple model and carried out stochastic 
particle simulations. Here we consider a simple one-dimensional reaction-diffusion system with three chemicals 
(Xi, X2, and X 3 ) and four reactions: 

fci 

X2 + X3 — ► X2 + X\ , 
X 3 + X\ — ^ 2X3, 

2Xi — X\ + Xj. 

Here, we assume fci, &2 ^> &3 > fc.4. We take k\ = k% = br, k 3 = ar, and = r (r > 0, 1 < a -C b) for further 
discussion. We assume that the system is closed with regards to the molecules X m (m — 1, 2, 3). Thus, N, the 
total number of molecules (or c, the total concentration), is conserved. 

In the continuum limit, each Ci(t, x), the concentration of Xi at time t and position x, obeys the following 
reaction-diffusion equation: 

9ci , , * , o 2\ t-. 9 2 C! , . 

— = -br(c 1 -C2)c 3 -r(c 1 -ac 2 ) + D 1 -^-, (7) 

— = r( Cl -ac 2 ) + D 2 -^, (8) 

— = Mci - C2 ) C3+A5 — (9) 

where is the diffusion constant of Xi. For simplicity, we assume Di = D for all i. This reaction-diffusion 
equation has homogeneous fixed point solutions with (ci, C2, C3) = (0, 0, c), {\facj (-y/ci+ 1), c/ (\/a + 1), 0) for all 
2;. By linear stability analysis, it is straightforward to show that only the former is stable. Starting from any 
initial conditions, the partial-differential equation system simply converges to this stable fixed point. 

In this system, the chemical X\ is produced by X2 molecules. (In the notation of §2, i — 2 and j = 1.) If i\ 1 
the Kuramoto length of Xi, is shorter than the average distance between X2 molecules, X\ is localized around 
the X2 molecules, as discussed in §2. Then, the reaction 2X\ — > X\ + X 2 , which is at second order of Xi, is 
accelerated. Using eq. (0 in §2, we obtain the acceleration factor 

° = 1 + 27kfi (10) 

On the other hand, the lifetime of X 2 is so long that X2 is not localized. Thus, the reaction 2X 2 — + X2 + X\ is 
not accelerated. 

Now we study the self-consistent solution from c\ and C2, following the argument of §2. (By a path from 
X\ to X2, there is a direct feedback to i, i.e., p = 2(— i) in the notation of §2.) Here, we consider the case 
where iVi, N2 <C N3, so that C3 « c. Then, the average lifetime of X\ is about X/ibrc^) « l/(brc); we assume 
that l\ = \j2D/{brc) for further discussion. When N\, N 2 <C N 3 , the two reactions X 2 + X 3 — » X 2 + X\ 
and X 3 + X\ — > 2X3 are much faster than the other two and maintain N\ ~ N 2 . Then, the rates of the other 
reactions satisfy 

The rate of (2X 1 -> X x + X 2 ) _ ak 4 Nf _ a 

The rate of (2X 2 ^X 2 +X 1 )^ k 3 N 2 "a' ( ' 

When a = a, the two reactions are balanced, which leads to a novel fixed point. Assuming c\,C2 <C C3 and 
C3 = c, and following eq. H10|) . we obtain the condition for the balance 

ci=c 2 = X (=c). (12) 

2^n{a — l)£i 

Subsequently, we investigate the stability of this fixed point. For c\ = c s + 5c\ and C2 = c s + <5c2, we obtain 

a = 1 + (a ~ 1)Cs = a - — -5c 2 + 0(602) (13) 

C2 C s 
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Figure 3: (a) Time series of N% and (b) the average concentrations c\ and c 2 for several values of N. a = 4, b = 
100, D = 100, and L x = 1000. We fix rc = 1 (i.e., r = 1000/JV), so that l x = y/2D/{brc) = y/2. For N > 1000 
cases, N 2 converges to the stochastic fixed point. There, C2 should be l/(2-y/7r(a — l)^i) = 1/6v2tt « 0.066 (i.e., 
N 2 = c 2 L x ~ 66) for all cases, which agrees well with the simulation. When N = 500, the system reaches the 
unstable fixed point (2c/3, c/3, 0). With the reaction-diffusion equation, c\ and c 2 rapidly converge to (the 
dotted line "continuous" in (a) shows c 2 for c = 1). 



from cqs. <|10|) and l|12l) . We take into account the acceleration factor in the reaction-diffusion equation, we 
obtain 

(c'i \ / — 2ac s — be (3a — l)c s + be \ ( 5c\ \ 

62 ) \ 2ac s -(3a-l)c s J \ Sc 2 J 

+o(Sc 1 ,5o 2 ). (14) 

For any a, b > 1, this Jacobi matrix has two negative eigenvalues, implying that the fixed point is stable 1 . 

In the simulations, we have found that the system converges to the novel fixed point. In fact, we have 
measured c 2 at the fixed point for a certain a = a numerically. Figure0]shows the relation between l/c 2 t\ and 
a, which agrees rather well with the theoretical estimation in eq. (|10|l . 

In summary, we demonstrated numerically that the discreteness of molecules yields a novel stochastic steady 
state in a reaction-diffusion system, in agreement with the theoretical estimation. 



4 Coexistence of Domains with different Kuramoto Lengths 

In the example of the previous section, the spatial homogeneity is assumed at a coarse-grained level, and indeed, 
this homogeneous state was stable. However, due to fluctuations, some spatial inhomogeneity exists in SSS, and 
a domain that is deviated from SSS may be produced. Even if this deviated state is unstable in the continuum 
limit, it may be preserved over a long time, if the concentration of the molecule to destabilize it is so low that 
its discreteness is essential. If the average time to produce this deviated domain from SSS and the lifetime of 
the state is balanced, the two regions, SSS and the deviated state with different concentrations of molecules 
and Kuramoto lengths, may coexist. We give a simple example for it here. 

Again we consider the same reactions as the preceding section. Now, we assume that the diffusion of X3 
is slower than the others, and set (Di, D 2 , D3) = (100-D, 100D, D) (D > 0). In this model, there are two 
fixed-point states: one is the stochastic steady state mentioned above (which we call State A); the other is the 
unstable fixed point with C3 = (State B), besides the stable fixed point in the continuum limit (which does 
not appear here; see Tabled for the stability of the states mentioned here). In Fig. El we give an example of 
snapshot pattern of the model. In the figure, except several spots with large C3 that corresponds to the state A, 
most other regions fall onto the state B that should be unstable in the continuum limit. Indeed, this pattern 

1 It is also stable against spatially inhomogeneous perturbations for any D > 0. 
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Figure 4: The acceleration factor a = a, plotted against \jo%l\. b = 100, N — 1000, L x = 1000, sampled over 
5000 (r = 1) or 50000 (r = 0.1) time units. We measure the relation from average C2 in simulations. The result 
agrees with the theoretical estimation a = 1 + „ J: , well. 



state 


Ci 


continuous 


discrete 


the unstable fixed point of the R-D eq. 
the stochastic steady state 
the stable fixed point of the R-D eq. 


c 3 = 
ci,c 2 ,c 3 > 
c\ = c-i = 


unstable 
(not fixed point) 
stable 


unstable 

stable 
unstable 



Table 1 : Stability of each fixed point in the continuous reaction-diffusion equation and in the stochastic system 
taking into account discreteness in molecules. 

is not transient, and the fraction of the state B is stationary, in the long-term simulation. This suggests the 
possibility that the state B, unstable in the continuum limit, may be sustained over a finite period due to the 
discreteness in molecules, which forms a domain in space. The two states A and B coexist in space and form a 
domain structure. 

First, we consider stability of each of the states in more detail. The state B is unstable against the inflow 
of A3. If an A3 molecule enters into a region of state B (Region B), it can be amplified and form a new region 
(spot) of the state A (Region A). From linear stability analysis, we find that the degree of instability of the 
region B against the flow of A3 is proportional to (ci — Ca). Note that in the region A, c\ sa C2, while in the 
region B, c\ ~ \/~ac2, and C3 = 0. Here, the concentration ci is almost uniform in space because of its long 
lifetime. Thus, the degree of the instability of the Region A, that is the rate of growth of C3, depends mainly 
on the distribution of X\ (see Fig. EJ. 

Accordingly, the Kuramoto length of X\ is relevant to determine if the state B is invaded or not. Here 
it should be noted that in the region A, Ci is smaller than that in the region B. Hence in the vicinity of the 
region A within the domain of the state B, c\ is still small, as long as it is within the Kuramoto length of X\ 
of the region A. Thus the instability is weak there, which prevents a novel region A (A3 spot) growing in the 
vicinity of the existing region A. Hence, the interval between two neighboring regions A should be longer than 
the Kuramoto length of X\. Assuming that N 3 <C Ni, N 2 and (01,02,03) = (-y=^y, ^ +1 , 0) (i.e., the unstable 
fixed point of the reaction-diffusion equation) in the region B, we obtain the Kuramoto length of X\ in the 
region B as 

t / on 71 , w / 20Qg(V^+l) 

tlB = V 2£) l( fc 4Cl) —\ 7= • 

y r^Jac 

Since A3 can be amplified by using X\ in the region B, penetration of A3 molecules into the region B must 
be rare in order to maintain the region B. The penetration length is given by the Kuramoto length of A3, that 
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Figure 5: The distribution of each chemical (a snapshot at t — 1000). r — 1.28, a = 4, b = 100, D = 1, 
N = 3200, and L x = 4096. (For these parameters, £\b ~ 17 and £35 « 0.24.) Plotted distribution is obtained 
by averaging the molecule snapshot pattern, with the bin size Ax = 16. There appear spots of A3, in which 
the stochastic steady state with a > 1 is realized. There are very few A3 molecules between the spots, where 
the system stays around the unstable fixed point (2c/3, c/3, 0). 

is computed as 



for the region B. For 1 < a <C (1006) 2 , £33 <C £\b-, which implies that A3 molecules seldom reach the area where 
A3 is strongly amplified. Thus, the border of regions A and B is maintained for long time 2 . 

On the other hand, due to the fluctuation inherent in SSS, the molecule X3 may be extinct within some 
area of the region A, with some probability. Hence, the regions A and B coexist in space, as shown in Fig. El 
As shown, the region A is localized only as spots, and other parts are covered by the region B. 

Note that in the corresponding reaction-diffusion equations, the state A, the stochastic steady state, cannot 
be realized, while the state B is unstable. Indeed, the reaction-diffusion equation system is quickly homogenized 
and converges to the stable fixed point with c\ — C2 = 0. Hence both the regions A and B, as well as a domain 
structure from the two, can exist only as a result of the discreteness of molecules, and are immediately destroyed 
in the continuum limit. 

Note that in the present model, SSS does not have so-called the Turing instability 3 , and there is no char- 
acteristic wavelength. Still the spatial structure here has some characteristic length, as given by the minimal 
size of the region B, estimated by the Kuramoto length lis- In Fig. we have plotted the spatial power 
spectrum of the concentrations of A3 . Although there is no clear peak, there is a very broad increase around 
the wavenumber of 0.01, that corresponds to average domain size of the region B. 

In general, the discreteness of molecules can induce novel states not seen in the continuum limit. For example, 
in a randomly connected catalytic reaction network, there often exist several fixed points with some chemicals 
going extinct, when the number of molecules is small, while there is only one attractor in the continuum limit. 
These discreteness-induced states may coexist in space in a similar way as discussed above. Kuramoto length 
will be a useful index to determine the behavior around the border of the regions. 

5 Summary and Discussion 

In the present paper, we have reported a novel steady state in a system with reaction and diffusion, induced by 
discreteness in molecules. This state cannot be represented by a continuum description, i.e., partial differential 

2 For this reason, we set D3 relatively small. If D3 is larger, the border is blurred and the two regions are mixed. 
3 It is also possible that the acceleration of reactions by the discreteness induces or enhances the Turing instability in certain 
systems. 
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Figure 6: Schematic diagram of the border between the regions A and B. 
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Figure 7: Power spectra of the distribution of X 3 . r = 1, a = 4, b = 100, N = 16384, and L x = 16384. For 
small D, there appears a broad peak at the wave length around the order of 10 2 , which corresponds to the 
interval of the spots. For large D, X3 spreads over the space, and the peak disappears. 
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equation (reaction-diffusion equation), but is sustained by amplification of some reaction due to localization 
of some molecule Xj. This localization is possible if the molecule species that produce Xj is "discrete", in 
the sense that its average number within the Kuramoto length of the Xj molecule is less than 1. We have 
formulated a theory to obtain a self-consistent solution for the concentrations of Xi and Xj, in relationship 
with the amplification rate of the reaction involved with Xj . For some reaction system, there is a solution with 
amplification rate larger than 1, that leads to the existence of stochastic steady state due to the discreteness in 
molecule number. The stability of this solution is also computed within this theoretical formulation. 

We have also numerically studied a simple reaction-diffusion system to demonstrate the validity of the theory. 
Indeed, a novel stochastic steady state is observed, as predicted theoretically. We have also extended our theory 
to include the self-consistent determination of the Kuramoto length. Following this extension, we have provided 
a numerical example, to show formation of domains with different Kuramoto lengths. 

The alteration of the steady state by the localization, as well as our formulation for it is quite general. 
Provided that the conditions 

(i) Chemical Xi generates another chemical species Xj. 

(ii) The lifetime of Xj is short or the diffusion of Xj is slow so that the Kuramoto length of Xj is much 
smaller than the distance between Xi molecules. 

(iii) The localization of the molecule Xj accelerates some reactions. 

Then, the discreteness can alter the dynamics, from that by the continuum description. The last condition is 
satisfied if the second or higher order reaction is involved in the species Xj . Finally, if 

(iv) the acceleration of the reaction in (iii) alters the density of Xi molecules, through some reaction(s), 

the density of Xi is determined self-consistently with the acceleration factor, resulting in a novel steady state. 

Note that the localization effect by the discreteness of catalytic molecules itself is also noted by Shnerb et 
al. JOI- In their study, however, the density of the catalyst is fixed as an externally given value. Thus the 
concentration of the product, localized around the catalyst, diverges in time. In our theory, the density of the 
catalyst (Xi) changes autonomously and reaches a suitable value by following the discreteness effect. 

The self-consistent solution scheme to obtain this discreteness-induced stochastic state can be extended to 
a case with several components ...,ik, and the corresponding set of chemicals 31,32, ■■■,jk satisfying the 

conditions (i)-(iii). In such case, the feedback process in (iv) is not necessarily direct from j m to i m . If there is 
a feedback from the set of chemicals 31,32, ■ ■■,jk to the set i\, i%, %k (condition (iv)'), the above scheme for the 
self-consistent dynamics we presented here works. With this extension, there is a variety of possibilities, that 
can lead to stochastic rhythm or pattern formation induced by discreteness of molecules, which is not seen in 
the continuum limit. For example, in a catalytic reaction network with many components and a limited number 
of total molecules, there always exist several species that are minority in number, and the conditions (i)-(ii) are 
naturally satisfied, while with higher order catalytic reaction the condition (iii) is often satisfied. In this case, 
minority molecules become a key factor to determine a macroscopic state with rhythm or pattern. (Note in this 
case, other molecules can be abundant in number, or indeed it is better to have such abundant species, so that 
the stochastic state is stabilized.) 

In fact, biochemical reaction networks involve a huge number of species, while the total number of molecules 
is not necessarily so large. In a cell, lots of chemicals work at low concentration in the order of 1 nM or less. The 
diffusion is sometimes restricted, surrounded by macro-molecules, and may be slow. In such an environment, it 
is probable that the average distance between the molecules of a given chemical species is much larger than the 
Kuramoto lengths of some other chemical species. Some chemicals are localized around some other molecules. 
Furthermore, biochemical systems contain various higher order reactions (for example, catalyzed by enzyme 
complexes). In conjunction with the localization, such reactions can be accelerated. Hence the conditions (i)- 
(iii) are ubiquitously satisfied in intra-cellular biochemical reaction networks. In addition, since the biochemical 
reactions involve complex feedback process through mutual catalytic networks, the condition (iv) or (iv)' is 
naturally satisfied. 

Accordingly, it will be important to study the amplification of some reaction and its maintenance through 
feedback will be relevant to biochemical reactions. Indeed, some molecules that are minority in number some- 
times play a key role in biological function. Relevance of minority molecules is also discussed from the viewpoint 
on a control mechanism of a cell, in relationship with the kinetic origin of information 1X3 • 
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The importance of our theory is not restricted to biological problems. Verification of our result will be possible 
by suitably designing a reaction system, with the use of, say, microreactors or vesicles. The acceleration and 
maintenance of some reactions by localization of molecules will be important to design some function in such 
micro-reactor systems. 
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